
#include "Covariance.hpp"
#include "nr.h"

Covariance::Covariance ( const NRMat<double>& p_EigenVec ) : 
  m_EigenVec(p_EigenVec), m_EigenVal(0.0,p_EigenVec.nrows()), 
  m_OffDiag(0.0,p_EigenVec.nrows()), m_SortVec(p_EigenVec.nrows())
{
  int n=m_SortVec.size();
  for ( int i=0; i<n; ++i )
    m_SortVec[i] = i;
}
Covariance::~Covariance() { }

void Covariance::PCA()
{
  // Householder transformation
  NR::tred2 (m_EigenVec, m_EigenVal, m_OffDiag);
  // QL algorithm
  // EigenVec becomes the eigenvector matrix
  // EigenVal becomes the eigenvalues
  NR::tqli(m_EigenVal, m_OffDiag, m_EigenVec);
  //copy eigenvalues & eigenvectors
  //NRVec<double> EigenVal(m_EigenVal);
  //NRMat<double> EigenVec(m_EigenVec);
  // Find the top eigenvalues
  // Sort Eigenvalues
  //NR::eigsrt(m_SortVec,EigenVal,EigenVec);
  EigenSort(m_SortVec,m_EigenVal);
}

void Covariance::EigenSort(NRVec<int> &s,const NRVec<double> &ev)
{
  int i,j,k,sp;
  double p;
  NRVec<double> d(ev);

  int n=d.size();
  for (i=0; i<n-1; i++) {
    p=d[k=i];
    for ( j=i; j<n; j++ )
      if (d[j] >= p)
	p=d[k=j];
    if (k != i) {
      d[k]=d[i];
      d[i]=p;
      sp=s[k]; // keep track of changes
      s[k]=s[i];
      s[i]=sp;
    }
  }
}
